AD-4037  318 


UNCLASSIFIED 


NAVAL  HEALTH  RESEARCH  CENTtR  SAN  DIEGO  CALIF  F/S  14/9 

A SET  OF  POP-12  PROGRAMS  FOR  COMPUTING  SPECTRA  OF  POWER • COHERE— ETC im 
Dec  71  CM  NUTE 

73-4  ml 


i °m 

AA0373I8 


DATE 

FILMED 

4^77 


^ SET  OF  PDp-12  ^ROGRAMS  FOR  ^OMPUTING  JPECTRA 
OF  POWER,  COHERENCE,  AND  PHASE  . 


H./NUTE 


REPORT  NO/  73-6 


/JUVAL  HEALTH 


SAN  DIEGO,  CALIFORNIA  92152 


NAVAL  MEDICAL  RESEARCH  AND  DEVELOPMENT  COMMAND 

BETH6SOA,  MARYLAND 

£?/  C</3l 


A SET  OF  POP-12  PROGRAMS  FOR  COMPUTING 
SPECTRA  OF  POUER,  COHERENCE,  AND  PHASE 


Cyril  H.  flute 

Navy  Medical  Neuropsychiatric  Research  Unit 
San  Diego,  California 


ABSTRACT 

i A time-varying  function  may  be  analyzed  for  Its  frequency  con- 
tent by  means  of  either  the  simple  Fourier  transformation  or 
the  smoothed  power  density  spectrum.  Two  or  more  functions 
may  be  related  through  coherence  and  phase  spectra.  Working 
from  either  analog  voltages  or  previously  prepared  digital 
tapes,  a set  of  programs  is  available  for  making  these  compu- 
tations. Runs  have  been  made  with  artificially  generated 
pseudo-random  numbers,  and  the  results  are  presented  to  indi- 
cate the  sampling  variability  to  be  expected  under  certain 
conditions. 


INTRODUCTION 
Deterministic  Functions 

Two  types  of  frequency  analysis  may  be  used  with  a 
sequence  of  digital  measurements  from  a time-varying 
function.  The  first  is  classical  Fourier  analysis. 
For  example,  assume  an  epoch  length  of  16  seconds 
with  a sampling  rate  of  128  Hz.  This  provides  a set 
of  2048  raw  input  data  points,  from  which  a Fourier 
transformation  will  obtain  1024  sine  coefficients 
and  1024  cosine  coefficients.  These  may  be  con- 
verted to  the  amplitudes  and  phases  of  a set  of  pure 
sinusoids  ranging  from  zero  to  64  Hz  by  increments 
of  .0625,  (1/16)  Hz.  The  basic  formula  states  that 
the  minimum  frequency  and  the  sharpest  frequency 
resolution  that  are  available  are  equal  to  the 
reciprocal  of  the  epoch  length.  The  maximum  fre- 
quency that  can  be  studied  is  half  the  sampling  fre- 
quency, known  as  the  Nyquist  or  folding  frequency. 

For  a deterministic  waveform,  like  the  AC  voltage  on 
a domestic  power  line,  such  an  analysis  might  be 
useful  for  measuring  the  line  frequency,  since  any 
other  sample  Is  likely  to  give  results  that  are  vir- 
tually Identical.  - For  discriminations  finer  than 
.0625  Hz,  there  Is  no  alternative  but  to  take  longerj 
samples.  With  a fixed  computer  buffer  size,  a 
trade-off  must  be  made  between  the  range  of  frequen- 
cies to  be  studied  and  the  maximum  frequency  resolu- . 
tion. 

Random  Functions 

An  entirely  different  class  of  functions  consists  of[ 
stationary  random,  or  stochastic  processes.  By  def-; 
inltlon,  these  do  not  Include  any  pure  sinusoids  but' 
may  Include  energy  that  tends  to  lie  within  certain  ' 
'limited  frequency  Intervals.  For  example,  a human 
being  with  a high  alpha  wave  will  have  a number  of 
strong  Fourier  components  In  the  neighborhood  of  10 
Hz,  but  no  two  runs  will  ever  yield  an  Identical  set 
of  amplitudes.  To  obtain  meaningful  (l.e.  consis-  i 
tent  and  unbiased)  results,  one  must  resort  to  some  I 
form  of  statistical  averaging.  The  squared  ampli- 
tudes, proportional  to  either  numerical  variance  or 
electrical  power,  of  all  of  the  frequencies  within 
successive  frequency  ranges  are  lumped  together  to 
give  relatively  gross  estimates  of  the  magnitudes  1 
and  locations  of  the  dominant  components.  The  I 


frequencies  originally  obtained  by  simple  Fourier 
transformation  have  been  called  the  “elementary,"  or 
"line"  frequencies,  due  to  the  appearance  of  the 
discrete  plots  of  their  squared  amplitudes.  If  a 
number,  1C,  of  these  values  are  added  together,  the 
sum  has  a chi-squared  distribution  with  2K  degrees 
of  freedom.1  This  function  Is  the  third  trade-off 
to  be  considered  in  setting  up  a spectral  analysis. 
If  T is  the  epoch  length  in  seconds,  and  S the  digi - 
tlzTng  rate  in  Hz,  the  following  relations  hold: 

(1)  Elementary  frequency  = 

(2)  Maximum  frequency  * 

(3)  No.  of  measurements  = 

(4)  Spectral  resolution  ■ 

(5)  Degrees  of  freedom  • 

In  brief,  the  analysis  may  be  designed  to  maximize 
the  confidence  associated  with  either  the  height  or 
the  position  of  a power  peak  along  the  frequency 
axis,  but  both  goals  cannot  be  achieved  simultane- 
ously. This  basic  uncertainty  has  been  shown  to  be 
mathematically  equivalent  to  the  Heisenberg  princi- 
ple In  atomic  physics.2  Similar  uncertainties  apply 
to  estimates  of  coherence  and  phase.  These  are 
derived  from  the  complex  Wlshart  distribution,  which 
is  an  extension  of  the  chi -squared. 3’1' 

, PROGRAM  OPERATION 

! 

The  routines  described  here,  which  have  also  been 
submitted  to  DECUS,  are  built  around  the  Fast 
Fourier  Transform  for  Real  Variables,  written  by 
1 James  Rothman  and  distributed  under  the  DECUS  number 
,8-143.  Originally  written  for  a PDP-8  without 
[Extended  Arithmetic  Element,  It  was  titled  FFT-R. 

[The  version  now  used  Is  FFT-E,  still  entirely  In  P- 
mode,  but  using  the  additional  hardware  and  with  a > 
minor  round-off  error  corrected  In  the  multiplica- 
tion subroutine.  The  coding  added  by  the  present 
[author  has  been  arranged  to  avoid  the  need  for  any 
[changes  to  the  Rothman  package  as  distributed. 

Two  versions  have  been  written  to  accept  either  ana- 
log voltages  «Ma  the  A/0  converter,  or  digital 
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measurements  via  LINCtape.  Because  of  the  needs  of 
our  laboratory,  the  tape  Input  version  Is  the  one 
that  has  been  more  fully  tested  and  submitted  to 
DECUS.  The  Internal  parameters  are  set  up  for  16 
seconds  of  data,  sampled  at  128  Hz,  and  requiring  an 
output  frequency  resolution  of  ,5  Hz.  However, 
these  are  easily  converted,  and  the  results  of  doing 
so  in  one  instance  are  described  below.  In  any 
case,  2048  points  from  one  channel  are  stored  In  the 
upper  4K  of  core,  while  the  FFT  routine  operates  on 
an  equal  number  of  points  from  the  other  channel  In 
lower  core.  The  routine  overlays  the  Input  data 
with  its  computed  Fourier  coefficients. 

After  one  channel  has  been  analyzed,  the  contents  of 
the  two  buffers  are  Interchanged,  and  the  second 
channel  Is  analyzed.  From  the  two  sets  of  Fourier 
coefficients,  the  two  auto-power  spectra  and  the 
real  and  imaginary  parts  of  the  cross-power  spectrum 
are  computed  and  stored  in  double  precision  Integer 
form  on  an  intermediate  -output  UNCtape  (Figure  1). 
All  of  these  computations  take  about  20-25  seconds. 

In  the  analog  version.  Intended  for  on-line  opera- 
tion, the  KW12  clock  Is  running  with  an  Interrupt 
routine  that  digitizes  Jhe  voltage  jnputs^  While  a 
30-second  batch  of  data  is  being  analyzed,  a suc- 
ceeding epoch  Is  being  collected  and  temporarily 
stored  on  LIHCtape  on  Unit  0.  As  soon  as  the  inter- 
mediate results  have  been  written  out  on  Unit  1,  the 
temporary  file  on  Unit  0 Is  unloaded  and  the  cycle 
restarted.  The  process  Is  entirely  automatic,  with 
successive  groups  of  one,  two,  or  three  blocks  writ- 
ten on  Unit  1 with  each  iteration,  depending  on  the 
parameters  inserted  to  the  program. 

The  01AL  program  that  does  this  work  Is  filed  under 
the  name  XS  for  the  tape  input  version  or  XSPECTRM 
for  the  analog  version.  After  a set  of  runs  has 
been  completed,  a F0CAL-1 2 program,  $XS,**  is  called 
In  and  executed  for  each  epoch  to  be  analyzed. 
Because  of  storage  limitations  on  an  8K  machine, 
overlays  have  to  be  used,  with  programs  $XS  and 
JTYPEXS  automatically  overlaying  each  other.  They 
produce  a final  typed  listing  of  frequency,  two 
auto-power  spectra,  and  coherence  and  phase. 

I 

•FOCAl-12  programs  are  distinguished  by  the  Initial 
$ In  their  names. 


SAMPLE  VARIABILITY 

Unfortunately,  spectrum  analysis  is  not  a simple 
tool  that  can  be  used  without  a good  understanding 
of  the  underlying  statistical  principles.  Although 
the  theoretical  distributions  of  computed  sample 
values  have  been  worked  out  on  the  basis  of  certain 
convenient  mathematical  assumptions,1 *3  one  cannot 
be  quite  certain  of  their  applicability  to  real 
problems.  Besides,  since  coherence  and  phase  are 
jointly  distributed  according  to  rules  that  are  con- 
siderably less  familiar  to  most  users  than  the  chi- 
squared  distribution,  special  care  must  be  exercised 
lln  their  use.  Some  simulations  have  been  made  with 
Gaussian  pseudo-random  number  sequences  Input  as 
data  In  order  to  determine  the  sampling  distribu- 
tions of  the  computed  parameters. 

I* 

Figure  2 shows  the  program  output  when  there  Is  zero 
correlation  between  the  two  input  channels.  With  an 
.epoch  length  of  16  .seconds  and  an  output  frequency 
resolution  of  .5  Hz,  each  of  the  estimated  power 
densities,  coherences,  and  phases  has  16  degrees  of 
freedom.  Because  of  some  shortcuts  taken  in  setting 
up  these  runs,  the  total  variances  printed  at  the 
top  of  the  sheet  are  not  in  the  same  units  as  the 
power  densities  given  below.  However,  Instead  of 
being  uniformly  distributed,  as  expected,  they  are 
seen  to  range  all  the  way  from  approximately  200 
units  all  the  way  to  approximately  1200  units.  The 
coherences,  which  would  all  be  zero  If  they  were 
consistent  and  unbiased  estimators  of  squared  corre- 
lation as  a function  of  frequency,  are  seen  to  go  as 
high  as  0.183.  Now  If  the  same  data  are  run  with  a 
broader  frequency  resolution  that  has  64  degrees  of 
freedom  (see  Figure  3),  the  random  distribution  of 
estimators  covers  a much  more  narrow  range.  The 
power  densities  vary  by  no  more  than  a factor  of 
two,  and  the  maximum  reported  coherence  is  five 
times  less  than  In  the  previous  run. 

Two  more  runs  were  made  under  similar  conditions, 
except  that  the  simulated  input  data  sequences  had  a 
known  squared  cross-correlation  of  0.50  (Figures  4 
and  5).  The  distribution  of  the  estimated  power 
densities  Is  about  the  same  as  before,  as  one  might 
expect,  but  the  computed  coherences  are  much  higher. 
In  Figure  2,  with  only  16  degrees  of  freedom,  the 
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Figure  1:  Steps  In  spectral  and  cross-spectral  analysis 


minimum  coherence  is  less  than  the  reported  maximum 
where  there  was  no  correlation  between  the  channels. 

It  appears  that  with  no  more  than  16  degrees  of 
freedom,  one  cannot  with  complete  certainty  distin- 
guish populations  with  zero  coherence  from  ones  with 
0.50  coherence.  However,  there  is  still  information 
to  be  gained  from  such  a computer  run,  because  In 
the  0.50  case,  most  of  the  computed  coherences  are 
well  above  the  one  reported  in  the  zero  case.  With 
64  degrees  of  freedom,  of  course,  the  distinction 
between  the  two  populations  is  immediate  and  obvious. 

CONCLUSION 

I. 

A set  of  programs  for  doing  spectral  and  cross- 
spectral  analyses  has  been  submitted  to  DECUS.  They 
are  designed  for  input  from  LINCtape  and  assume  a 
sampling  rate  of  128  Hz  over  a period  of  16  seconds, 
with  an  output  frequency  resolution  of  0.50  Hz. 
However,  other  versions  exist  for  analog  input  in 
real  time  with  different  frequency  and  epoch  param- 
eters. The  author  will  be  glad  to  advise  In  adapt- 
ing the  submitted  programs  to  the  needs  of  a partic- 
ular user.  These  programs  were  run  with  simulated 
data  and  have  shown  that  with  16  degrees  of  freedom, 
one  can  distinguish  fairly  reliably  between  popula- 
tions whose  true  coherence  is  zero  and  those  with  a 
true  coherence  of  0.50.  With  64  degrees,  the  dis- 
tinction Is  much  more  certain,  and  even  finer  reso- 
lutions are  possible. 
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Figure  4:  results  of  a run  similar  to  that  in  Fig.  2,  except  that  the 
squared  correlation  coefficient  between  the  two  input  series 
Number  of  degrees  of  freedom  ■ 16. 
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